*------------------------------------------------------------------------------*
*	Paper:			International Trade and Public Protest
*					Evidence from Russian Regions
*	Authors:		Tabea Palmtag, Tobias Rommel, Stefanie Walter
*	Version:		June 10, 2020
*	Task:			Empirical Analysis
*------------------------------------------------------------------------------*


clear all
version 14
set mat 3000
set more off
capture log close



***** 1. Dataset
*------------------------------------------------------------------------------*

*** Set directory
local dir = "/Users//`=c(username)'/Downloads/PRW_2020_ISQ_Replication/"
cd "`dir'"

*** Load data
use "PRW_2020_ISQ_regiondata.dta", clear

*** Graph preferences
* ssc install blindschemes
* set scheme plotplain



***** 2. Data Preparation
*------------------------------------------------------------------------------*

*** Convert regional GDP from RUB to USD
gen grp=reg_grp*(xr_rus)^(-1)
lab var grp "Gross regional product (in million USD)"

*** Trade openness
* Imports and exports
gen tradeopen=((reg_exporttofor+reg_exporttosng+ ///
	reg_importtofor+reg_importtosng)/grp)*100
lab var tradeopen "Imports+Exports (% of GDP)"
gen lntrade = ln(tradeopen + 1)
lab var lntrade "Imports+Exports (% of GDP), ln"
* Exports
gen exportopen=((reg_exporttofor+reg_exporttosng)/grp)*100
lab var exportopen "Exports (% of GDP)"
gen lnexport = ln(exportopen + 1)
lab var lnexport "Exports (% of GDP), ln"
* Imports
gen importopen=((reg_importtofor+reg_importtosng)/grp)*100
lab var importopen "Imports (% of GDP)"
gen lnimport = ln(importopen + 1)
lab var lnimport "Imports (% of GDP), ln"

*** FDI openness
gen fdiopen=((reg_fdi_incl_total/1000)/grp)*100
lab var fdiopen "FDI inflows (% of GDP)"
gen lnfdi = ln(fdiopen + 1)
lab var lnfdi "FDI inflows (% of GDP), ln"

*** Generate average education level (secondary + tertiary)
* Secondary and tertiary (main models)
bysort rcode: egen edshare=mean(higheduc) if year>=2007
gen educshare = edshare-34
* Tertiary (main models)
bysort rcode: egen edshareter=mean(tereduc) if year>=2007
gen educshareter = edshareter-17
* Secondary + 2*tertiary (main models)
bysort rcode: egen edshareind=mean(higheduc+tereduc) if year>=2007
gen educshareind = edshareind-54
* Secondary and tertiary (robustness checks)
bysort rcode: egen edshare2=mean(higheduc) if year>=2003
gen educshare2 = edshare2-36
* Secondary and tertiary (regional mechanism)
bysort rcode: egen edshare4=mean(higheduc) if year>=2000
gen educshare4 = edshare4-36

*** Generate GRP growth
sort rcode year
by rcode: gen reg_grpgr=(grp-grp[_n-1])/grp[_n-1]

*** Population size in millions
replace reg_pop =reg_pop/1000000

*** Population density
gen reg_density = reg_pop/reg_area

*** Logarithm of GRP per capita
gen lngrppc = ln(reg_grp_pc)

*** Newspaper coverage
gen lnnews = ln(reg_newspaper+1)

*** Road density
gen lnroadden = ln(reg_autoroadden)

*** Distance to Moscow
gen lndistmos = ln(reg_disttom+1)



***** 3. Descriptive Statistics
*------------------------------------------------------------------------------*

*** Sample size
xtset rcode year
xtnbreg protest_ikd2 cl.lntrade##c.educshare l.lnfdi l.reg_pop ///
	l.reg_urbanshare l.lngrppc l.reg_grpgr l.reg_levelofunempl ///
	l.lnnews l.lndistmos l.lnroadden, re
gen sample = e(sample)==1

*** Table A1 in the Online Appendix
tabstat protest_ikd2 protest_ikd_econ protest_mma protest_icews protest_kprf2 ///
	lntrade educshare lnfdi reg_pop reg_urbanshare lngrppc reg_grpgr ///
	reg_levelofunempl lnnews lndistmos lnroadden rents pressfill ///
	KPRFmandateshare pctRussian if sample==1, statistics(N mean sd min max) ///
	format(%9.2f)

*** Correlation between protest indicators
corr protest_ikd2 protest_ikd_econ protest_mma protest_icews protest_kprf2 ///
	if sample==1

*** Histograms
hist protest_ikd2 if sample==1, ///
	xtitle("Number of grassroots IKD protests") ///
	graphregion(fcolor(white) ilcolor(white) lcolor(white) color(white) ///
		ifcolor(white) style(none)) ///
	name(hist_ikd, replace) fysize(40)
hist tradeopen if sample==1, ///
	xtitle("Exposure to international trade") ///
	graphregion(fcolor(white) ilcolor(white) lcolor(white) color(white) ///
		ifcolor(white) style(none)) ///
	name(hist_trade, replace) fysize(40)

*** Prepare trade and FDI exposure
preserve
* Keep protest and trade variables
keep rname rid year protest_ikd2 protest_kprf2 tradeopen lntrade fdiopen lnfdi
keep if year>=2000 & year!=2014
* Generate merge variable for individual-level data
gen conversion = 	.		
replace conversion = 	1	if rname == 	"Altai K"
replace conversion = 	2	if rname == 	"Amur"
replace conversion = 	3	if rname == 	"Chelyabinsk"
replace conversion = 	4	if rname == 	"Chuvash"
replace conversion = 	5	if rname == 	"Kabardino-Balkar"
replace conversion = 	6	if rname == 	"Kaliningrad"
replace conversion = 	7	if rname == 	"Kaluga"
replace conversion = 	8	if rname == 	"Komi"
replace conversion = 	9	if rname == 	"Krasnodar"
replace conversion = 	10	if rname == 	"Krasnoyarsk"
replace conversion = 	11	if rname == 	"Kurgan"
replace conversion = 	12	if rname == 	"Leningrad"
replace conversion = 	13	if rname == 	"Lipetsk"
replace conversion = 	14	if rname == 	"Moscow C"
replace conversion = 	15	if rname == 	"Moscow O"
replace conversion = 	16	if rname == 	"Nizhny Novgorod"
replace conversion = 	17	if rname == 	"Novosibirsk"
replace conversion = 	18	if rname == 	"Orenburg"
replace conversion = 	19	if rname == 	"Penza"
replace conversion = 	20	if rname == 	"Perm"
replace conversion = 	21	if rname == 	"Primorsky"
replace conversion = 	22	if rname == 	"Rostov"
replace conversion = 	23	if rname == 	"Saint Petersburg"
replace conversion = 	24	if rname == 	"Saratov"
replace conversion = 	25	if rname == 	"Smolensk"
replace conversion = 	26	if rname == 	"Stavropol"
replace conversion = 	27	if rname == 	"Tambov"
replace conversion = 	28	if rname == 	"Tatarstan"
replace conversion = 	29	if rname == 	"Tomsk"
replace conversion = 	30	if rname == 	"Tula"
replace conversion = 	31	if rname == 	"Tyumen"
replace conversion = 	32	if rname == 	"Udmurt"
replace conversion = 	33	if rname == 	"Volgograd"
save "glob_stats.dta", replace
* Prepare data for maps
keep if year>=2007 & year<=2012
collapse (first) rname (sum) protest_ikd2 protest_kprf2 (mean) tradeopen ///
	lntrade fdiopen lnfdi, by(rid)
save "map_stats.dta", replace

*** Maps
* ssc install shp2dta
* ssc install geo2xy
* ssc install spmap
* Load shapefile
shp2dta using "PRW_2020_ISQ_RUSadm1.shp", ///
	database("map_data") ///
	coordinates("map_coord") ///
	genid(id) replace
* Prepare coordinates file
use "map_coord", clear
replace _X = 180 if _X<0
geo2xy _Y _X if _X>=-180 & _X<=180, replace projection(albers)
drop if _Y<=4000000
save "map_coord", replace
* Prepare region identifiers
use "map_data", clear
gen rname = NAME_1
replace rname = 	"Adygea"	if rname == 	"Adygey"
replace rname = 	"Altai K"	if rname == 	"Altay"
replace rname = 	"Altai R"	if rname == 	"Gorno-Altay"
replace rname = 	"Arkhangelsk"	if rname == 	"Arkhangel'sk"
replace rname = 	"Buryatia"	if rname == 	"Buryat"
replace rname = 	"Chechen"	if rname == 	"Chechnya"
replace rname = 	"Chukotka"	if rname == 	"Chukot"
replace rname = 	"Ingushetia"	if rname == 	"Ingush"
replace rname = 	"Kabardino-Balkar"	if rname == 	"Kabardin-Balkar"
replace rname = 	"Kalmykia"	if rname == 	"Kalmyk"
replace rname = 	"Khakassia"	if rname == 	"Khakass"
replace rname = 	"Khanty-Mansi"	if rname == 	"Khanty-Mansiy"
replace rname = 	"Mari El"	if rname == 	"Mariy-El"
replace rname = 	"Magadan"	if rname == 	"Maga Buryatdan"
replace rname = 	"Moscow C"	if rname == 	"Moscow City"
replace rname = 	"Moscow O"	if rname == 	"Moskva"
replace rname = 	"Nizhny Novgorod"	if rname == 	"Nizhegorod"
replace rname = 	"North Ossetia-Alania"	if rname == 	"North Ossetia"
replace rname = 	"Oryol"	if rname == 	"Orel"
replace rname = 	"Primorsky"	if rname == 	"Primor'ye"
replace rname = 	"Saint Petersburg"	if rname == 	"City of St. Petersburg"
replace rname = 	"Ulyanovsk"	if rname == 	"Ul'yanovsk"
replace rname = 	"Yamalo-Nenets"	if rname == 	"Yamal-Nenets"
replace rname = 	"Zabaykalsky"	if rname == 	"Zabaykal'ye"
replace rname = 	"Astrakhan"	if rname == 	"Astrakhan'"
replace rname = 	"Perm"	if rname == 	"Perm'"
replace rname = 	"Ryazan"	if rname == 	"Ryazan'"
replace rname = 	"Stavropol"	if rname == 	"Stavropol'"
replace rname = 	"Tver"	if rname == 	"Tver'"
replace rname = 	"Tyumen"	if rname == 	"Tyumen'"
replace rname = 	"Yaroslavl"	if rname == 	"Yaroslavl'"
replace rname = 	"Jewish Autonomous"	if rname == 	"Yevrey"
drop if rname==		"Sverdlovsk" & CCA_1!="0"
save "map_data", replace
* Merge descriptives
use "map_stats", clear
merge 1:1 rname using "map_data"
drop _merge
spmap protest_ikd2 using "map_coord" if rname!="Kaliningrad", ///
	id(id) clnumber(9) fcolor(Greys2) ///
	ocolor(white white white white white white white white white) ///
	osize(vthin vthin vthin vthin vthin vthin vthin vthin vthin) ///
	legenda(off) title("Grassroots IKD Protest") name(map_ikd, replace)
spmap tradeopen using "map_coord" if rname!="Kaliningrad", ///
	id(id) clnumber(9) fcolor(Greys2) ///
	ocolor(white white white white white white white white white) ///
	osize(vthin vthin vthin vthin vthin vthin vthin vthin vthin) ///
	legenda(off) title("International Trade Exposure") name(map_trade, replace)

*** Figure 2 in the Main Text
graph combine map_ikd map_trade hist_ikd hist_trade, col(2) ysize(3) xsize(5.5)
graph save "figure_descriptives", replace
restore


	
***** 4. Protest models
*------------------------------------------------------------------------------*

*** Model specification for IKD and KPRF
* Dataset structure
xtset rcode year
* Baseline specifications
global X1 cl.lntrade  c.educshare
global X2 cl.lntrade##c.educshare
global C1 l.reg_pop l.reg_urbanshare l.lngrppc l.reg_grpgr ///
	l.reg_levelofunempl l.lnnews l.lndistmos l.lnroadden reg_density l.lnfdi 
global C2 rents pressfill KPRFmandateshare pctRussian 
* Sample size
xtnbreg protest_ikd2 $X2 $C1, re
gen sample1 = e(sample)==1

*** Model specification for MMA and ICEWS
* Model setup
global X3 cl.lntrade  c.educshare2
global X4 cl.lntrade##c.educshare2
* Sample size
xtnbreg protest_mma $X4 $C1, re
gen sample2 = e(sample)==1

*** Regression models
* IKD protest, unconditional
xtnbreg protest_ikd2 $X1 $C1 if sample1==1, re
est store ikd1
* IKD protest, conditional
xtnbreg protest_ikd2 $X2 $C1 if sample1==1, re
est store ikd2
* IKD protest, conditional + extended controls
xtnbreg protest_ikd2 $X2 $C1 $C2 if sample1==1, re
est store ikd3
* IKD protest, conditional + extended controls + lagged protest
xtnbreg protest_ikd2 $X2 $C1 $C2 l.protest_ikd2 if sample1==1, re
est store ikd4
* IKD protest, conditional without onshore oil
gen oilregion = 0
replace oilregion = 1 if rname=="Adygea" | rname=="Astrakhan" | ///
	rname=="Bashkortostan" | rname=="Chukotka" | rname=="Irkutsk" | ///
	rname=="Kaliningrad" | rname=="Kalmykia" | rname=="Kamchatka" | ///
	rname=="Khabarovsk" | rname=="Khanty-Mansi" | rname=="Kirov" | ///
	rname=="Komi" | rname=="Krasnodar" | rname=="Krasnoyarsk" | ///
	rname=="Novosibirsk" | rname=="Omsk" | rname=="Rostov" | ///
	rname=="Sakha" | rname=="Sakhalin" | rname=="Saratov" | ///
	rname=="Sverdlovsk" | rname=="Tatarstan" | rname=="Tyumen" | ///
	rname=="Volgograd" | rname=="Yamalo-Nenets"
xtnbreg protest_ikd2 $X2 $C1 if sample1==1 & oilregion==0, re
est store ikd5
* IKD economic protest, conditional
xtnbreg protest_ikd_econ $X2 $C1 if sample1==1, re
est store ikd_econ
* MMA protest, conditional
xtnbreg protest_mma $X2 $C1 if sample2==1, re
est store mma
* ICEWS protest, conditional
xtnbreg protest_icews $X2 $C1 if sample2==1, re
est store icews
* KPRF protest, conditional
xtnbreg protest_kprf2 $X2 $C1 if sample1==1, re
est store kprf
* KPRF economic protest, conditional
xtnbreg protest_kprf_econ2 $X2 $C1 if sample1==1, re
est store kprf_econ
* Table 1 in Main Text
esttab ikd1 ikd2 ikd3 ikd4 ikd5 ikd_econ mma icews kprf, ///
	cells(b(star fmt(%9.3f)) se(par fmt(2))) style(fixed) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label ///
	stats(N N_g p ll aic, labels("# of observations" "# of regions" ///
		"Prob > Chi2" "Log likelihood" "AIC") fmt(0 0 3 2 2)) ///
	order(*lntrade* educ*)

*** Additional control variables
* IKD protest, including FDI interaction
xtnbreg protest_ikd2 $X2 $C1 cl.lnfdi##c.educshare  if sample1==1, re
est store control1
* IKD protest, including cost of living
xtnbreg protest_ikd2 $X2 $C1 reg_costliving if sample1==1, re
est store control2
* IKD protest, including % of people living below subsistence
xtnbreg protest_ikd2 $X2 $C1 reg_belowcost if sample1==1, re
est store control3
* IKD protest, including cost of a fixed basket of goods
xtnbreg protest_ikd2 $X2 $C1 reg_gdsfixed if sample1==1, re
est store control4
* IKD protest, including economic inequality
xtnbreg protest_ikd2 $X2 $C1 reg_mincgini if sample1==1, re
est store control5
* IKD protest, including interregional migration
xtnbreg protest_ikd2 $X2 $C1 reg_interregmigration if sample1==1, re
est store control6
* IKD protest, including infant mortality
xtnbreg protest_ikd2 $X2 $C1 reg_infantmort if sample1==1, re
est store control7
* IKD protest, including life expectancy
xtnbreg protest_ikd2 $X2 $C1 reg_lifeexp if sample1==1, re
est store control8
* IKD protest, including public expenditure differential index
xtnbreg protest_ikd2 $X2 $C1 reg_ibr if sample1==1, re
est store control9
* Table A4 in Online Appendix
esttab control1 control2 control3 control4 control5 control6 control7 control8 control9, ///
	cells(b(star fmt(%9.3f)) se(par fmt(2))) style(fixed) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label ///
	stats(N N_g p ll aic, labels("# of observations" "# of regions" ///
		"Prob > Chi2" "Log likelihood" "AIC") fmt(0 0 3 2 2)) ///
	order(*lntrade* educ*)

*** Additional robustness checks
* Longer average education
xtnbreg protest_ikd2 $X4 $C1 $C2 if sample1==1, re
est store robust1
* Tertiary education only
xtnbreg protest_ikd2 cl.lntrade##c.educshareter $C1 $C2 if sample1==1, re
est store robust2
* Secondary + 2*tertiary education
xtnbreg protest_ikd2 cl.lntrade##c.educshareind $C1 $C2 if sample1==1, re
est store robust3
* Exports only
xtnbreg protest_ikd2 cl.lnexport##c.educshare $C1 $C2 if sample1==1, re
est store robust4
* Imports only
xtnbreg protest_ikd2 cl.lnimport##c.educshare $C1 $C2 if sample1==1, re
est store robust5
* Fixed effects negative binomial (without secondary education share)
xtnbreg protest_ikd2 cl.lntrade cl.lntrade#c.educshare $C1 if sample1==1, fe
est store robust6
* Federal district fixed effects
nbreg protest_ikd2 $X2 $C1 $C2 i.dname if sample1==1
est store robust7
* Economic region fixed effects
nbreg protest_ikd2 $X2 $C1 $C2 i.ename if sample1==1
est store robust8
* Small regions only (< median)
gen smallregionm = 0
replace smallregion = 1 if reg_area<=75
xtnbreg protest_ikd2 $X2 $C1 $C2 if sample1==1 & smallregion==1, re
est store robust9
* Table A5 in Online Appendix
esttab robust1 robust2 robust3 robust4 robust5 robust6 robust7 robust8 robust9, ///
	cells(b(star fmt(%9.3f)) se(par fmt(2))) style(fixed) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label ///
	stats(N N_g p ll aic, labels("# of observations" "# of regions" ///
		"Prob > Chi2" "Log likelihood" "AIC") fmt(0 0 3 2 2)) ///
	order(*lntrade* educ*)

*** Figure 3 in Main Text
xtnbreg protest_ikd2 $X2 $C1 $C2 if sample1==1, re
est store m1
margins, dydx(l.lntrade) at(educshare=(0(1)30))
marginsplot, level(90) ///
	plotopts(msymbol(i) lcolor(black)) ///
	ciopts(recast(rline) lpattern(dash) lcolor(black)) ///
	title("") ///
	xtitle("Share of population with secondary education (in %)", ///
		size(large) margin(medsmall)) ///
	ytitle("Marginal effect of trade exposure", size(vlarge) margin(small)) ///
	ylabel(, angle(hori) labsize(medsmall) nogrid) ///
	xscale(noline range(-.5 26.5)) yscale(noline) ///
	xlabel(1 "35" 6 "40" 11 "45" 16 "50" 21 "55" 26 "60", ///
		labsize(medsmall) nogrid) ///
	yline(0, lcolor(black) lwidth(.2)) ///
	graphregion(fcolor(white) ilcolor(white) lcolor(white) color(white) ///
		ifcolor(white) style(none)) ///
	name(ikd_trade, replace) xsize(3) ysize(2) ///
	addplot(hist educshare, bin(21) yaxis(2) legend(off) ///
		ytitle("", axis(2)) ylabel(none, axis(2)) ///
		xlabel(1 "35" 6 "40" 11 "45" 16 "50" 21 "55" 26 "60", ///
			labsize(medsmall) nogrid) ///
		fcolor(none) lcolor(gs12))
graph save "figure_effect", replace



***** 5. Mechanism: regional level
*------------------------------------------------------------------------------*

*** Model specification
xtset rcode year
global X cl.lntrade##c.educshare4 
global C l.reg_pop l.reg_urbanshare l.lngrppc l.reg_grpgr ///
	l.reg_levelofunempl l.lnnews l.lndistmos l.reg_density l.lnroadden

*** Average wage level (positive interaction)
gen lnwage = ln(reg_avwage/reg_grp)
xtpcse lnwage $X $C if year>=2000, pair
est store mech1

*** Consumption per capita (positive interaction)
gen cons = PERSCONSUMPC/1000
xtpcse cons $X $C if year>=2000, pair
est store mech2

*** Level of employment (positive interaction)
gen empl = (EMPL/(reg_pop*1000))
xtpcse empl $X $C if year>=2000, pair
est store mech3

*** Table 2 in Main Text
esttab mech1 mech2 mech3, ///
	cells(b(star fmt(%9.3f)) se(par fmt(2))) style(fixed) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label ///
	stats(N N_g r2 p, labels("# of observations" "# of regions" ///
		"R squared " "Prob > Chi2") fmt(0 0 2 3)) ///
	order(*lntrade* educ*)



***** 6. Mechanism: individual risk
*------------------------------------------------------------------------------*

*** Russian survey panel
* Load data
use "PRW_2020_ISQ_RLMSdata.dta", clear
* Restrict to 2000-2013
drop if year<2000
drop if year>2013
* Keep adults only
drop if child==1

*** Merge trade and FDI exposure
* Generate merge variable
gen conversion = 	.		
replace conversion = 	1	if region == 	58
replace conversion = 	1	if region == 	84
replace conversion = 	2	if region == 	93
replace conversion = 	3	if region == 	106
replace conversion = 	3	if region == 	107
replace conversion = 	4	if region == 	48
replace conversion = 	5	if region == 	77
replace conversion = 	6	if region == 	67
replace conversion = 	7	if region == 	14
replace conversion = 	8	if region == 	89
replace conversion = 	8	if region == 	105
replace conversion = 	9	if region == 	9
replace conversion = 	9	if region == 	129
replace conversion = 	10	if region == 	66
replace conversion = 	10	if region == 	73
replace conversion = 	11	if region == 	46
replace conversion = 	12	if region == 	1
replace conversion = 	13	if region == 	72
replace conversion = 	14	if region == 	138
replace conversion = 	14	if region == 	140
replace conversion = 	15	if region == 	142
replace conversion = 	16	if region == 	116
replace conversion = 	17	if region == 	161
replace conversion = 	18	if region == 	47
replace conversion = 	19	if region == 	117
replace conversion = 	20	if region == 	12
replace conversion = 	21	if region == 	92
replace conversion = 	22	if region == 	137
replace conversion = 	23	if region == 	141
replace conversion = 	24	if region == 	70
replace conversion = 	24	if region == 	100
replace conversion = 	25	if region == 	135
replace conversion = 	26	if region == 	52
replace conversion = 	27	if region == 	33
replace conversion = 	28	if region == 	45
replace conversion = 	29	if region == 	71
replace conversion = 	30	if region == 	136
replace conversion = 	31	if region == 	86
replace conversion = 	32	if region == 	10
replace conversion = 	33	if region == 	39
label define conversion ///					
	1 "Altai K"	2 "Amur" 3 "Chelyabinsk" 4 "Chuvash" 5 "Kabardino-Balkar" ///
	6 "Kaliningrad"	7 "Kaluga" 8 "Komi" 9 "Krasnodar" 10 "Krasnoyarsk" ///
	11 "Kurgan"	12 "Leningrad" 13 "Lipetsk"	14 "Moscow C" 15 "Moscow O" ///
	16 "Nizhny Novgorod" 17 "Novosibirsk" 18 "Orenburg" 19 "Penza" 20 "Perm" ///
	21 "Primorsky" 22 "Rostov" 23 "Saint Petersburg" 24 "Saratov" ///
	25 "Smolensk" 26 "Stavropol" 27 "Tambov" 28 "Tatarstan" 29 "Tomsk" ///
	30 "Tula" 31 "Tyumen" 32 "Udmurt" 33 "Volgograd"		
label values conversion conversion					

*** Merge trade and FDI
merge m:m conversion year using "glob_stats.dta"
drop if _merge==2
drop _merge

*** Education level
gen educlvl = diplom-1 if diplom<=6

*** Dependent variables
* Satisfaction with economic condition
gen satecon = j66_1
replace satecon = . if j66_1>5
* Concern about getting necessities
gen necess = j66
replace necess = . if j66>5
* Rank on social ladder
gen econlad = j62
replace econlad = . if j62>9
gen powlad = j63
replace powlad = . if j63>9
gen ladder = (econlad + powlad)/2
gen ladder_rev = (ladder*(-1))+10
* Unemployed and would like to find work
gen unemp = 0 if j1<=5
replace unemp = 1 if j1==5
gen unempwork = 0 if unemp!=.
replace unempwork = 1 if unemp==1 & j81==1

*** Control variables
* Gender
gen gender = 0 if h5==1
replace gender = 1 if h5==2
* Age
replace age = . if age>104
* Income
gen income = j60
replace income = . if j60>2265000
gen lnincome = ln(income+1)
* Marriage
gen married = 0 if j72_17<=6 | j322<=6
replace married = 1 if j72_17==2 | j322==2 | j322==3
* Second job
gen secjob = j32 if j32<=2
replace secjob = 0 if secjob==2
replace secjob = 0 if unemp==1
* Works for government
gen govfirm = j23 if j23<=2
replace govfirm = 0 if govfirm==2
replace govfirm = 0 if unemp==1
* Works for foreign firm
gen forfirm = j24 if j24<=2
replace forfirm = 0 if forfirm==2
replace forfirm = 0 if unemp==1
* Works for Russian firm
gen rusfirm = j25 if j25<=2
replace rusfirm = 0 if rusfirm==2
replace rusfirm = 0 if unemp==1
* Self-employed
gen selfemp = j26 if j26<=2
replace selfemp = 0 if selfemp==2
replace selfemp = 0 if unemp==1
* Pension
gen pension = j73 if j73<=2
replace pension = 0 if pension==2

*** Model specification
global X1 c.educlvl##c.lntrade i.year
global X2 c.educlvl##c.lnfdi i.year
global C1 gender age lnincome married i.urban secjob selfemp govfirm
xtset idind year

*** Regression models
* Economic dissatisfaction
meoprobit satecon $X1 $C1 if pension!=1 || idind:
est store satecon1
meoprobit satecon $X1 lnfdi $C1 if pension!=1 || idind:
est store satecon2
* Concerned about getting necessities
gen necess_rev = (necess*(-1))+6
meoprobit necess_rev $X1 $C1 if pension!=1 || idind:
est store necess1
meoprobit necess_rev $X1 lnfdi $C1 if pension!=1 || idind:
est store necess2
* Unemployed and wants to find work
meprobit unempwork $X1 $C1 if pension!=1 || idind:
est store unempwork1
meprobit unempwork $X1 lnfdi $C1 if pension!=1 || idind:
est store unempwork2
* Rank on social ladder
mixed ladder_rev $X1 $C1 if pension!=1 || idind:
est store ladder1
mixed ladder_rev $X1 lnfdi $C1 if pension!=1 || idind:
est store ladder2

*** Table A2 in the Online Appendix
tabstat unempwork necess_rev satecon ladder_rev lntrade educlvl lnfdi ///
	gender age lnincome married secjob selfemp govfirm urban ///
	if pension!=1, statistics(N mean sd min max) format(%9.2f)

*** Table 3 in Main Text
esttab unempwork2 necess2 satecon2 ladder2, ///
	cells(b(star fmt(%9.3f)) se(par fmt(2))) style(fixed) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label ///
	stats(N N_g p ll aic, labels("# of observations" "# of regions" ///
		"Prob > Chi2" "Log likelihood" "AIC") fmt(0 0 3 2 2)) ///
	order(*lntrade* educ*)

*** Figure
meprobit unempwork $X1 lnfdi $C1 if pension!=1 || idind:
margins, dydx(lntrade) at(educlvl=(0(1)5))
marginsplot, level(95) ///
	plotopts(msymbol(i) lcolor(black)) ///
	ciopts(recast(rline) lpattern(dash) lcolor(black)) ///
	title("") ///
	xtitle("Individual education level", size(large) margin(medsmall)) ///
	ytitle("Marginal effect of trade exposure", size(vlarge) margin(small)) ///
	ylabel(, angle(hori) labsize(medsmall) nogrid) ///
	xscale(noline range(-.5 5.5)) yscale(noline) ///
	xlabel(0 `" "Primary" "or less" "' 1 `" "Primary" "completed" "' ///
		2 `" "Lower" "secondary" "' 3 `" "Higher" "secondary" "' ///
		4 `" "Secondary" "vocational" "' 5 `" "Tertiary" "completed" "', ///
			labsize(medsmall) nogrid) ///
	yline(0, lcolor(black) lwidth(.2)) ///
	graphregion(fcolor(white) ilcolor(white) lcolor(white) color(white) ///
		ifcolor(white) style(none)) ///
	name(unemp_trade, replace) ///
	addplot(hist educlvl, discrete yaxis(2) legend(off) ///
		ytitle("", axis(2)) ylabel(none, axis(2)) ///
		xlabel(0 `" "Primary" "or less" "' 1 `" "Primary" "completed" "' ///
			2 `" "Lower" "secondary" "' 3 `" "Higher" "secondary" "' ///
			4 `" "Secondary" "vocational" "' 5 `" "Tertiary" "completed" "', ///
				labsize(medsmall) nogrid) ///
		fcolor(none) lcolor(gs12))

*** Figure
meoprobit necess_rev $X1 lnfdi $C1 if pension!=1 || idind:
margins, dydx(lntrade) at(educlvl=(0(1)5)) predict(outcome(5))
marginsplot, level(95) ///
	plotopts(msymbol(i) lcolor(black)) ///
	ciopts(recast(rline) lpattern(dash) lcolor(black)) ///
	title("") ///
	xtitle("Individual education level", size(large) margin(medsmall)) ///
	ytitle("Marginal effect of trade exposure", size(vlarge) margin(small)) ///
	ylabel(, angle(hori) labsize(medsmall) nogrid) ///
	xscale(noline range(-.5 5.5)) yscale(noline) ///
	xlabel(0 `" "Primary" "or less" "' 1 `" "Primary" "completed" "' ///
		2 `" "Lower" "secondary" "' 3 `" "Higher" "secondary" "' ///
		4 `" "Secondary" "vocational" "' 5 `" "Tertiary" "completed" "', ///
			labsize(medsmall) nogrid) ///
	yline(0, lcolor(black) lwidth(.2)) ///
	graphregion(fcolor(white) ilcolor(white) lcolor(white) color(white) ///
		ifcolor(white) style(none)) ///
	name(necess_trade, replace) ///
	addplot(hist educlvl, discrete yaxis(2) legend(off) ///
		ytitle("", axis(2)) ylabel(none, axis(2)) ///
		xlabel(0 `" "Primary" "or less" "' 1 `" "Primary" "completed" "' ///
			2 `" "Lower" "secondary" "' 3 `" "Higher" "secondary" "' ///
			4 `" "Secondary" "vocational" "' 5 `" "Tertiary" "completed" "', ///
				labsize(medsmall) nogrid) ///
		fcolor(none) lcolor(gs12))
	
*** Figure
meoprobit satecon $X1 lnfdi $C1 if pension!=1 || idind:
margins, dydx(lntrade) at(educlvl=(0(1)5)) predict(outcome(5))
marginsplot, level(95) ///
	plotopts(msymbol(i) lcolor(black)) ///
	ciopts(recast(rline) lpattern(dash) lcolor(black)) ///
	title("") ///
	xtitle("Individual education level", size(large) margin(medsmall)) ///
	ytitle("Marginal effect of trade exposure", size(vlarge) margin(small)) ///
	ylabel(, angle(hori) labsize(medsmall) nogrid) ///
	xscale(noline range(-.5 5.5)) yscale(noline) ///
	xlabel(0 `" "Primary" "or less" "' 1 `" "Primary" "completed" "' ///
		2 `" "Lower" "secondary" "' 3 `" "Higher" "secondary" "' ///
		4 `" "Secondary" "vocational" "' 5 `" "Tertiary" "completed" "', ///
			labsize(medsmall) nogrid) ///
	yline(0, lcolor(black) lwidth(.2)) ///
	graphregion(fcolor(white) ilcolor(white) lcolor(white) color(white) ///
		ifcolor(white) style(none)) ///
	name(satecon_trade, replace) ///
	addplot(hist educlvl, discrete yaxis(2) legend(off) ///
		ytitle("", axis(2)) ylabel(none, axis(2)) ///
		xlabel(0 `" "Primary" "or less" "' 1 `" "Primary" "completed" "' ///
			2 `" "Lower" "secondary" "' 3 `" "Higher" "secondary" "' ///
			4 `" "Secondary" "vocational" "' 5 `" "Tertiary" "completed" "', ///
				labsize(medsmall) nogrid) ///
		fcolor(none) lcolor(gs12))

*** Figure
mixed ladder_rev $X1 lnfdi $C1 if pension!=1 || idind:
margins, dydx(lntrade) at(educlvl=(0(1)5))
marginsplot, level(95) ///
	plotopts(msymbol(i) lcolor(black)) ///
	ciopts(recast(rline) lpattern(dash) lcolor(black)) ///
	title("") ///
	xtitle("Individual education level", size(large) margin(medsmall)) ///
	ytitle("Marginal effect of trade exposure", size(vlarge) margin(small)) ///
	ylabel(, angle(hori) labsize(medsmall) nogrid) ///
	xscale(noline range(-.5 5.5)) yscale(noline) ///
	xlabel(0 `" "Primary" "or less" "' 1 `" "Primary" "completed" "' ///
		2 `" "Lower" "secondary" "' 3 `" "Higher" "secondary" "' ///
		4 `" "Secondary" "vocational" "' 5 `" "Tertiary" "completed" "', ///
			labsize(medsmall) nogrid) ///
	yline(0, lcolor(black) lwidth(.2)) ///
	graphregion(fcolor(white) ilcolor(white) lcolor(white) color(white) ///
		ifcolor(white) style(none)) ///
	name(ladder_trade, replace) ///
	addplot(hist educlvl, discrete yaxis(2) legend(off) ///
		ytitle("", axis(2)) ylabel(none, axis(2)) ///
		xlabel(0 `" "Primary" "or less" "' 1 `" "Primary" "completed" "' ///
			2 `" "Lower" "secondary" "' 3 `" "Higher" "secondary" "' ///
			4 `" "Secondary" "vocational" "' 5 `" "Tertiary" "completed" "', ///
				labsize(medsmall) nogrid) ///
		fcolor(none) lcolor(gs12))

*** Figure 4 in Main Text 
graph combine unemp_trade necess_trade satecon_trade ladder_trade, xcommon
graph save "figure_individual", replace

*** Folder maintenance
erase "glob_stats.dta"
erase "map_coord.dta"
erase "map_data.dta"
erase "map_stats.dta"



***** 7. Mechanism: individual grievances
*------------------------------------------------------------------------------*

*** World Values Survey (Russia)
* Load data
use "PRW_2020_ISQ_WVSdata.dta", clear

*** Protest participation
gen protestrec = 0 if E221B==2
replace protestrec = 1 if E221B==1

*** Grievance indicators
* Life satisfaction
gen lifesat = A170 if A170>0
replace lifesat = (lifesat*(-1))+11
* Financial satisfaction
gen econsat = C006 if C006>0
replace econsat = (econsat*(-1))+11
* Income inequality
gen inequal = E035 if E035>0
replace inequal = (inequal*(-1))+11
* Not enough food (wave 6 only)
gen nofood = H008_01 if H008_01>0
recode nofood (1=4) (2=3) (3=2) (4=1)
* No cash income (wave 6 only)
gen nocash = H008_04 if H008_04>0
recode nocash (1=4) (2=3) (3=2) (4=1)

*** Controls
* Education
gen educlvl = X025 if X025>0
recode educlvl (1=1) (2=2) (3=3) (5=3) (4=4) (6=4) (7=5) (8=6)
* Female
gen female = 0 if X001==1
replace female = 1 if X001==2
* Age
gen age = X003
* Married
gen married = 0 if X007==3 | X007==4 | X007==5 | X007==6
replace married = 1 if X007==1 | X007==2
* Employment status
gen unemp = X028==7
* Income deciles
gen income = X047 if X047>0
* Politics important
gen imppol = A004 if A004>0
recode imppol (1=4) (2=3) (3=2) (4=1)
* Federal district
gen region = X048WVS

*** Table A3 in the Online Appendix
tabstat protestrec lifesat econsat inequal nofood nocash female age income ///
	married unemp imppol educlvl, statistics(N mean sd min max) format(%9.2f)

*** Regression models
* Specification
global C female age income married unemp imppol ib6.educlvl wave6
* Models
xtprobit protestrec lifesat $C, i(region) vce(cluster region)
est store prot1
xtprobit protestrec econsat $C, i(region) vce(cluster region)
est store prot2
xtprobit protestrec inequal $C, i(region) vce(cluster region)
est store prot3
xtprobit protestrec nofood $C, i(region) vce(cluster region)
est store prot4
xtprobit protestrec nocash $C, i(region) vce(cluster region)
est store prot5
* Table 4 in Main Text
esttab prot1 prot2 prot3 prot4 prot5, ///
	cells(b(star fmt(%9.3f)) se(par fmt(2))) style(fixed) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label ///
	stats(N N_g r2_p p, labels("# of observations" "# of regions" ///
		"Pseudo R squared " "Prob > Chi2") fmt(0 0 2 3)) ///
	order(lifesat econsat inequal nofood nocash)

*** Regression models
* Specification
global C female age income married unemp imppol ib6.educlvl wave6 i.region
* Models
probit protestrec lifesat $C, cluster(region)
est store prot6
probit protestrec econsat $C, cluster(region)
est store prot7
probit protestrec inequal $C, cluster(region)
est store prot8
probit protestrec nofood $C, cluster(region)
est store prot9
probit protestrec nocash $C, cluster(region)
est store prot10
* Table
esttab prot6 prot7 prot8 prot9 prot10, ///
	cells(b(star fmt(%9.3f)) se(par fmt(2))) style(fixed) ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label ///
	stats(N N_g r2_p p, labels("# of individuals" "# of regions" ///
		"Pseudo R squared " "Prob > Chi2") fmt(0 0 2 3)) ///
	order(lifesat econsat inequal nofood nocash)



*------------------------------------------------------------------------------*


capture log close
exit
